home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / include / gsl / gsl_odeiv.h < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-20  |  9.4 KB  |  283 lines
  1. /* ode-initval/gsl_odeiv.h
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman
  21.  */
  22. #ifndef __GSL_ODEIV_H__
  23. #define __GSL_ODEIV_H__
  24.  
  25. #include <stdio.h>
  26. #include <stdlib.h>
  27.  
  28. #undef __BEGIN_DECLS
  29. #undef __END_DECLS
  30. #ifdef __cplusplus
  31. # define __BEGIN_DECLS extern "C" {
  32. # define __END_DECLS }
  33. #else
  34. # define __BEGIN_DECLS /* empty */
  35. # define __END_DECLS /* empty */
  36. #endif
  37.  
  38. __BEGIN_DECLS
  39.  
  40.  
  41. /* Description of a system of ODEs.
  42.  *
  43.  * y' = f(t,y) = dydt(t, y)
  44.  *
  45.  * The system is specified by giving the right-hand-side
  46.  * of the equation and possibly a jacobian function.
  47.  *
  48.  * Some methods require the jacobian function, which calculates
  49.  * the matrix dfdy and the vector dfdt. The matrix dfdy conforms
  50.  * to the GSL standard, being a continuous range of floating point
  51.  * values, in row-order.
  52.  *
  53.  * As with GSL function objects, user-supplied parameter
  54.  * data is also present. 
  55.  */
  56.  
  57. typedef struct  
  58. {
  59.   int (* function) (double t, const double y[], double dydt[], void * params);
  60.   int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);
  61.   size_t dimension;
  62.   void * params;
  63. }
  64. gsl_odeiv_system;
  65.  
  66. #define GSL_ODEIV_FN_EVAL(S,t,y,f)  (*((S)->function))(t,y,f,(S)->params)
  67. #define GSL_ODEIV_JA_EVAL(S,t,y,dfdy,dfdt)  (*((S)->jacobian))(t,y,dfdy,dfdt,(S)->params)
  68.  
  69.  
  70. /* General stepper object.
  71.  *
  72.  * Opaque object for stepping an ODE system from t to t+h.
  73.  * In general the object has some state which facilitates
  74.  * iterating the stepping operation.
  75.  */
  76.  
  77. typedef struct 
  78. {
  79.   const char * name;
  80.   int can_use_dydt_in;
  81.   int gives_exact_dydt_out;
  82.   void * (*alloc) (size_t dim);
  83.   int  (*apply)  (void * state, size_t dim, double t, double h, double y[], double yerr[], const double dydt_in[], double dydt_out[], const gsl_odeiv_system * dydt);
  84.   int  (*reset) (void * state, size_t dim);
  85.   unsigned int  (*order) (void * state);
  86.   void (*free)  (void * state);
  87. }
  88. gsl_odeiv_step_type;
  89.  
  90. typedef struct {
  91.   const gsl_odeiv_step_type * type;
  92.   size_t dimension;
  93.   void * state;
  94. }
  95. gsl_odeiv_step;
  96.  
  97.  
  98. /* Available stepper types.
  99.  *
  100.  * rk2    : embedded 2nd(3rd) Runge-Kutta
  101.  * rk4    : 4th order (classical) Runge-Kutta
  102.  * rkck   : embedded 4th(5th) Runge-Kutta, Cash-Karp
  103.  * rk8pd  : embedded 8th(9th) Runge-Kutta, Prince-Dormand
  104.  * rk2imp : implicit 2nd order Runge-Kutta at Gaussian points
  105.  * rk4imp : implicit 4th order Runge-Kutta at Gaussian points
  106.  * gear1  : M=1 implicit Gear method
  107.  * gear2  : M=2 implicit Gear method
  108.  */
  109.  
  110. #ifdef GSL_EXPORTS
  111. __declspec(dllexport) const gsl_odeiv_step_type *gsl_odeiv_step_rk2;
  112. #elif defined(GSL_IMPORTS)
  113. __declspec(dllimport) const gsl_odeiv_step_type *gsl_odeiv_step_rk2;
  114. #else
  115. extern const gsl_odeiv_step_type *gsl_odeiv_step_rk2;
  116. #endif
  117. #ifdef GSL_EXPORTS
  118. __declspec(dllexport) const gsl_odeiv_step_type *gsl_odeiv_step_rk4;
  119. #elif defined(GSL_IMPORTS)
  120. __declspec(dllimport) const gsl_odeiv_step_type *gsl_odeiv_step_rk4;
  121. #else
  122. extern const gsl_odeiv_step_type *gsl_odeiv_step_rk4;
  123. #endif
  124. #ifdef GSL_EXPORTS
  125. __declspec(dllexport) const gsl_odeiv_step_type *gsl_odeiv_step_rkf45;
  126. #elif defined(GSL_IMPORTS)
  127. __declspec(dllimport) const gsl_odeiv_step_type *gsl_odeiv_step_rkf45;
  128. #else
  129. extern const gsl_odeiv_step_type *gsl_odeiv_step_rkf45;
  130. #endif
  131. #ifdef GSL_EXPORTS
  132. __declspec(dllexport) const gsl_odeiv_step_type *gsl_odeiv_step_rkck;
  133. #elif defined(GSL_IMPORTS)
  134. __declspec(dllimport) const gsl_odeiv_step_type *gsl_odeiv_step_rkck;
  135. #else
  136. extern const gsl_odeiv_step_type *gsl_odeiv_step_rkck;
  137. #endif
  138. #ifdef GSL_EXPORTS
  139. __declspec(dllexport) const gsl_odeiv_step_type *gsl_odeiv_step_rk8pd;
  140. #elif defined(GSL_IMPORTS)
  141. __declspec(dllimport) const gsl_odeiv_step_type *gsl_odeiv_step_rk8pd;
  142. #else
  143. extern const gsl_odeiv_step_type *gsl_odeiv_step_rk8pd;
  144. #endif
  145. #ifdef GSL_EXPORTS
  146. __declspec(dllexport) const gsl_odeiv_step_type *gsl_odeiv_step_rk2imp;
  147. #elif defined(GSL_IMPORTS)
  148. __declspec(dllimport) const gsl_odeiv_step_type *gsl_odeiv_step_rk2imp;
  149. #else
  150. extern const gsl_odeiv_step_type *gsl_odeiv_step_rk2imp;
  151. #endif
  152. #ifdef GSL_EXPORTS
  153. __declspec(dllexport) const gsl_odeiv_step_type *gsl_odeiv_step_rk4imp;
  154. #elif defined(GSL_IMPORTS)
  155. __declspec(dllimport) const gsl_odeiv_step_type *gsl_odeiv_step_rk4imp;
  156. #else
  157. extern const gsl_odeiv_step_type *gsl_odeiv_step_rk4imp;
  158. #endif
  159. #ifdef GSL_EXPORTS
  160. __declspec(dllexport) const gsl_odeiv_step_type *gsl_odeiv_step_bsimp;
  161. #elif defined(GSL_IMPORTS)
  162. __declspec(dllimport) const gsl_odeiv_step_type *gsl_odeiv_step_bsimp;
  163. #else
  164. extern const gsl_odeiv_step_type *gsl_odeiv_step_bsimp;
  165. #endif
  166. #ifdef GSL_EXPORTS
  167. __declspec(dllexport) const gsl_odeiv_step_type *gsl_odeiv_step_gear1;
  168. #elif defined(GSL_IMPORTS)
  169. __declspec(dllimport) const gsl_odeiv_step_type *gsl_odeiv_step_gear1;
  170. #else
  171. extern const gsl_odeiv_step_type *gsl_odeiv_step_gear1;
  172. #endif
  173. #ifdef GSL_EXPORTS
  174. __declspec(dllexport) const gsl_odeiv_step_type *gsl_odeiv_step_gear2;
  175. #elif defined(GSL_IMPORTS)
  176. __declspec(dllimport) const gsl_odeiv_step_type *gsl_odeiv_step_gear2;
  177. #else
  178. extern const gsl_odeiv_step_type *gsl_odeiv_step_gear2;
  179. #endif
  180.  
  181.  
  182. /* Constructor for specialized stepper objects.
  183.  */
  184. gsl_odeiv_step * gsl_odeiv_step_alloc(const gsl_odeiv_step_type * T, size_t dim);
  185. int  gsl_odeiv_step_reset(gsl_odeiv_step * s);
  186. void gsl_odeiv_step_free(gsl_odeiv_step * s);
  187.  
  188. /* General stepper object methods.
  189.  */
  190. const char * gsl_odeiv_step_name(const gsl_odeiv_step *);
  191. unsigned int gsl_odeiv_step_order(const gsl_odeiv_step * s);
  192.  
  193. int  gsl_odeiv_step_apply(gsl_odeiv_step *, double t, double h, double y[], double yerr[], const double dydt_in[], double dydt_out[], const gsl_odeiv_system * dydt);
  194.  
  195. /* General step size control object.
  196.  *
  197.  * The hadjust() method controls the adjustment of
  198.  * step size given the result of a step and the error.
  199.  * Valid hadjust() methods must return one of the codes below.
  200.  *
  201.  * The general data can be used by specializations
  202.  * to store state and control their heuristics.
  203.  */
  204.  
  205. typedef struct 
  206. {
  207.   const char * name;
  208.   void * (*alloc) (void);
  209.   int  (*init) (void * state, double eps_abs, double eps_rel, double a_y, double a_dydt);
  210.   int  (*hadjust) (void * state, size_t dim, unsigned int ord, const double y[], const double yerr[], const double yp[], double * h);
  211.   void (*free) (void * state);
  212. }
  213. gsl_odeiv_control_type;
  214.  
  215. typedef struct 
  216. {
  217.   const gsl_odeiv_control_type * type;
  218.   void * state;
  219. }
  220. gsl_odeiv_control;
  221.  
  222. /* Possible return values for an hadjust() evolution method.
  223.  */
  224. #define GSL_ODEIV_HADJ_INC   1  /* step was increased */
  225. #define GSL_ODEIV_HADJ_NIL   0  /* step unchanged     */
  226. #define GSL_ODEIV_HADJ_DEC (-1) /* step decreased     */
  227.  
  228. gsl_odeiv_control * gsl_odeiv_control_alloc(const gsl_odeiv_control_type * T);
  229. int gsl_odeiv_control_init(gsl_odeiv_control * c, double eps_abs, double eps_rel, double a_y, double a_dydt);
  230. void gsl_odeiv_control_free(gsl_odeiv_control * c);
  231. int gsl_odeiv_control_hadjust (gsl_odeiv_control * c, gsl_odeiv_step * s, const double y0[], const double yerr[], const double dydt[], double * h);
  232. const char * gsl_odeiv_control_name(const gsl_odeiv_control * c);
  233.  
  234. /* Available control object constructors.
  235.  *
  236.  * The standard control object is a four parameter heuristic
  237.  * defined as follows:
  238.  *    D0 = eps_abs + eps_rel * (a_y |y| + a_dydt h |y'|)
  239.  *    D1 = |yerr|
  240.  *    q  = consistency order of method (q=4 for 4(5) embedded RK)
  241.  *    S  = safety factor (0.9 say)
  242.  *
  243.  *                      /  (D0/D1)^(1/(q+1))  D0 >= D1
  244.  *    h_NEW = S h_OLD * |
  245.  *                      \  (D0/D1)^(1/q)      D0 < D1
  246.  *
  247.  * This encompasses all the standard error scaling methods.
  248.  *
  249.  * The y method is the standard method with a_y=1, a_dydt=0.
  250.  * The yp method is the standard method with a_y=0, a_dydt=1.
  251.  */
  252.  
  253. gsl_odeiv_control * gsl_odeiv_control_standard_new(double eps_abs, double eps_rel, double a_y, double a_dydt);
  254. gsl_odeiv_control * gsl_odeiv_control_y_new(double eps_abs, double eps_rel);
  255. gsl_odeiv_control * gsl_odeiv_control_yp_new(double eps_abs, double eps_rel);
  256.  
  257.  
  258. /* General evolution object.
  259.  */
  260. typedef struct {
  261.   size_t dimension;
  262.   double * y0;
  263.   double * yerr;
  264.   double * dydt_in;
  265.   double * dydt_out;
  266.   double last_step;
  267.   unsigned long int count;
  268.   unsigned long int failed_steps;
  269. }
  270. gsl_odeiv_evolve;
  271.  
  272. /* Evolution object methods.
  273.  */
  274. gsl_odeiv_evolve * gsl_odeiv_evolve_alloc(size_t dim);
  275. int gsl_odeiv_evolve_apply(gsl_odeiv_evolve *, gsl_odeiv_control * con, gsl_odeiv_step * step, const gsl_odeiv_system * dydt, double * t, double t1, double * h, double y[]);
  276. int gsl_odeiv_evolve_reset(gsl_odeiv_evolve *);
  277. void gsl_odeiv_evolve_free(gsl_odeiv_evolve *);
  278.  
  279.  
  280. __END_DECLS
  281.  
  282. #endif /* __GSL_ODEIV_H__ */
  283.